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Abstract 

Using a cellular automata model, we simulate the British Govern- 
ment Policy (BGP) in the 2001 foot and mouth epidemic in Great 
Britain. When clinical symptoms of the disease appeared on a farm, 
there is mandatory slaughter (culling) of all livestock on an infected 
premise (IP). Those farms that neighbor an IP (contiguous premise, 
CP), are also culled, aka nearest neighbor interaction. Farms where 
the disease may be prevalent from animal, human, vehicle or airborne 
transmission (dangerous contact, DC), are additionally culled, aka 
next-to-nearest neighbor iteractions and lightning factor. The result- 
ing mathematical model possesses a phase transition, whereupon if 
the physical disease transmission kernel exceeds a critical value, catas- 
trophic loss of animals ensues. The non-local disease transport prob- 
ability can be as low as .01% per day and the disease can still be in 
the high mortality phase. We show that the fundamental equation 
for sustainable disease transport is the criticality equation for neutron 
fission cascade. Finally, we calculate that the percentage of culled 
animals that are actually healthy is ~ 30%. 
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1 Prologue 



On or about February 19th, 2001, on a swill feeding farm (the index farm) 
at Hcdon-on-the-Wall in the county of Tyne and Wear, type O Pan-Asian 
strain foot and mouth virus disease (FMD) was introduced into the United 
Kingdom. From there the virus spread (probably by wind) to a sheep farm 7 
km away, and then to 96 different locations in Great Britain. By the time the 
epidemic was identified, 43 farms were infected. To confront the epidemic, 
the British Government instituted a policy of culling (killing of animals) to 
eradicate the infection. One year later with over 4 million animals destroyed 
(with unquantifiable human misery). Great Britain was declared disease- free. 

2 Introduction 

Foot and mouth disease is an extremely contagious trans-species (Bovidae 
such as cattle, zebus, domestic buffaloes and yaks; sheep, swine, goats, all 
wild ruminants, and suidae) virus that can live temporally outside the host 
mammal. In this paper, we simulate in a mathematical model the 2001 BGP 
in controlling and finally eradicating foot and mouth virus in Great Britain. 
One could call the BGP "risk aversion" in that it relied on no vaccinations 
and culling of IP when clinical symptoms appeared. The main interest here 
is the question of whether the BGP risked an unacceptable loss of animals by 
not knowing that the underlying FMD mathematical model exists in different 
phases. 

The FMD stochastic disease propagation cellular automata model used 
in this paper has strong resemblance to the famous "forest fire" models^. In 
general, probabilistic cellular automata models are a modern tool for sim- 
ulations of dynamical systems. They can model complex spatio-temporal 
structures involving cooperative, time-dependent phenomena^. In contrast 
to the mentioned forest fire diffusion models which have forest regeneration, 
no steady-state is present with the BGP: the disease becomes eradicated 
(The British Government prohibited the restocking of culled farms - ALL 
livestock killed - until after the epidemic was officially declared over). Thus, 
the only question is whether the peculiar epidemiology of FMD causes the 
BGP to incur surprisingly severe losses. How this may happen is illustrated 
in Fig. 1. 
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Figure 1: Culling destroys the animals of farm A and adjoining farm B, but 
FMD can have already spread from farm B 

There are three time scales: incubation time, contagious time and clinical 
symptom time. All three time scales depend statistically on the host type of 
animal, as well as the particular strain of FMD. For purposes of illustration, 
let us choose a 3 day incubation time and a eight day clinical symptom time. 
Refering to Fig. 1, farm A is infected on day 1 and enters its contagious period 
on day 4. After day 3 onwards, it can transmit the disease to neighbors (farm 
B) with some probablity per day value, which depends on the type of farm 
B animals, their population and their distance away from farm A. On day 8, 
clinical symptoms appear at farm A, and according to BGP, it becomes an 
IP, with all animals culled (next day). Also, according to BGP, farm B is a 
CP, and so its animals are culled, as well. As can be seen from Fig. 1, there is 
a time period called the "time-lag" , whereby the animals at farm B may have 
been contagious and, before being culled, could have transmitted the disease 
to yet another farm. In this illustration, the possible time-lag is 3 days (an 
extreme case). Thus the culling of farms A and B does not necessarily stop 
the diffusion of the disease starting originally from farm A. The existence of 
the time-lag is the characteristic epidemiological feature of FMD. In the next 
section, examples of physical probabilities 11 are given from the FMD data, 
and show that the disease has a long spatial tail, making time-lags inevitable 
when a large number of infections are present. We will see that the cellular 
automata model of BGP with time-lag has a phase transition, or critical 
point, called the percolation limit. If the physical disease transmission kernel 
exceeds a critical value, the disease percolates throughout the landscape, in 
spite of the culling policy. 

Within this paper, three variants of the cellular automata model are pre- 
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sented, and discussed in more detail later on. All model variants use [1] a 
300 X 300 cellular landscape; [2] a conditional probability per day infection 
(n) which falls off with distance. This cellular automata fall-off distribution 
is shown in Fig. 2. 
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Figure 2: Cellular diffusion probability distribution 

The middle cell row = i, column = j is the IP. Those cells having a P 
diffusion probability (discussed below) are the CP. The cells having a P/3 
diffusion probability are the DC. The infection can jump from an IP to CP, 
DC or beyond. 

It is important to understand the difference between P and 11. When an 
IP is identified in the contagious state, other cells containing healthy animals 
are at risk from this IP. This probability risk increases with the number of 
healthy animals in the other cell and decreases with the distance from the 
IP. Thus n is the conditional probability per day: given that an IP exists 
X km away, what is the probability per day of infection 11 if the healthy 
farm contains N number of animals. The diffusion probability P, however, 
measures the probability per day that a cell will obtain the disease, and so 
it is equal to the probability that animals are actually present in the cell 
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(Panimais) times the probability per day that those animals will contact the 
disease (11): 

P Panimals X 11 (l) 

Since England has about 93,000 square miles of area, the 300 x 300 lattice 
has a physical mapping metric of about 1 square mile per cell. 

The first model variant, called the homogeneous model, has farms in 
all cells, with a uniform, conditional probability per day of infection from 
neighbors, U — constant. That is, an IP can infect nearest neighbors with 
a uniform, conditional probability H, if these neighbors still have animals. 
An IP also can infect second-nearest neighbor farms with a smaller, but still 
uniform, conditional probability 11/3, and third- nearest neighbors with a still 
smaller uniform, conditonal probability n/9. 

The second variant, called the heterogeneous model, adds a random spa- 
tial density of farms (labeled by the symbol u, where u can range from 
to 1). We use a farm density value u = 0.75. Thus, a cell has a 100 x i/ 
probability of having a farm. As with the homogeneous model, an IP can 
infect nearest neighbors with a uniform, conditional probability 11, if these 
neighbors still have animals, second-nearest neighbor farms with a uniform, 
conditional probability 11/3; and third-nearest neighbors with a uniform, 
conditonal probability 11/9. 

The third and most realistic BGP simulation has a heterogeneous spatial 
distribution of farms, a stochastic conditional probability per day infection 
of neighbors, 11 7^ constant, and, most importantly, a lightning factor (LF 
7^ 0). The LF represents the non-local transmission probability of FMD^. 



Model Variants 


n 


Farm density {i>) 


LF 


homogeneous case 


constant 


1 





heterogeneous case 


constant 


0.75 





BGP variant 


variable 


0.75 


constant (7^ 0) 



Table 1: The three models presented 
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3 Conditional Infection Probabilities per Day 



The extensive spread of the 2001 FMD in Great Britain allowed for the first 
time identification of the physical transmission probability kernel^. This ker- 
nel is an empirically derived function (depending on the susceptibility factor 
of each species of animal, the number of animals of each species and the 
distance from an IP) using FMD data taken directly from the 2001 epidemic. 
Using this kernel for cattle, we plot in Fig. 3 the infection probability per 
day for small farms (500 and 1000 cattle) and for large farms (2500 and 5000 
cattle). This probability per day of infection means that a farm containing 
N number of healthy animals situated X kms away from an IP, has an empir- 
ically derived probability per day of getting the disease. So taking the 5000 
cattle graph as an example: each day, a large farm containing 5000 healthy 
cattle has a 60% chance becoming infected with FMD if an IP is located 
6 km away. As already mentioned, wc use the symbol 11 to designate this 
conditional probability per day attribute associated with a farm containing 
healthy animals. 



The developed cellular automata model is a fully object oriented C-I--I- sim- 
ulation. The square lattice 300x300 = 90,000 cells (with row i and col- 
umn j) has both time and state information. The possible state information 
is: Cell[i]|j][l] = is culled (originally with animals, but now empty), 1 is 
diseased, 2 has unaffected animals, -1 means empty (no original animals). 
Cell[i][j][0] has the time information: On the day of the infection, Ccll[i][j][0] 
has the value: -(time-lag)- (incubation); for each time unit (a day) 1 is added 
to this such that when it is -(time-lag), contagious transmission can occur 
and when zero, culling about Ccll[i][j] occurs. An infected cell remains in- 
fected until culling occurs. Empty cells are just an impediment to the disease 
diffusion, not being capable of transmitting the disease. One of the more 
interesting questions that a simulation can address is to compute the per- 
centage of culled animals that arc actually healthy. It is relatively simple to 
keep track of the state of animals at the time of culling, and wc find that the 




4 Cellular Automata Model 
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Figure 3: Computed conditional probability per day of infection 11 vs. dis- 
tance (km), based on actual transmission kernel data^ from the Great Britain 
outbreak. 

percentage healthy can be as high as 40%, with a typical average of around 
30%. 

In Fig. 4, we give the logic of the program, while in Fig. 2, the individual 
diffusion probabilities per day (P) are shown, relative to IP Cell[i][j]. When 
IP Cell[i][j] has reached clinical symptoms, all CP are culled (those with dif- 
fusion probability P) and all DC are culled (those with diffusion probability 
P/3). 

Since the cellular automata model is a Monte Carlo simulation involving 
millions of random numbers, it is important to use an unbiased random 
number generator. In this paper, the "gold plated"^ ran2 was used. Each 
run had different random numbers because the C++ time( ) function was 
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For each timestep 
For each Cell[i][j] 

-> Check if lightning fector causes disease 

-> Check if cell is infected (IP) and ready to be culled 

> If true, cull nearest neighbor (CP) 

> If true, cull next-to-nearest neighbors (DC) 

-> Check if cell is infected and can transmit the virus 

> If true, roll probability for nearest neighbors 
with healthy animals 

> If true, roll probability for next-to-nearest 
neighbors with healthy animals 

-> Update each infected cell's time history 
-> Check if simulation can stop (no live virus in 
landscape) 

Figure 4: Program logic 

used to set the value of the random number seed. 

As mentioned previously, the main interest in the research is the question 
of whether the BGP risked an unacceptable loss of animals by not knowing 
that the underlying FMD mathematical model exists in different phases. 
Thus we are looking for critical behavior as we vary the model parameters. 
There are multiple ways to detect a phase transition. One is to plot the 
mortality or survivor ratio versus the transmission probability per day and 
see a sharp transition. A second method is to plot the duration time of 
the epidemic versus transmission probability per day and see that a peak 
duration appears. The larger the lattice size, the more pronounced is the 
phase transition. 

Once the existence of the phase transition is established, there are in- 
teresting quantities that can be investigated, such as critical exponents and 
the disease time correlation function. A particularly challenging problem is 
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to compute analytically the critical probability of the phase transition for 
variants of the model. In actual practice, one needs to go to extremely large 
lattices to see power law and fractal scaling behavior, because small lattices 
produce a finite width for the phase transition. The 300 x 300 lattice used 
here, though producing a substantial width to the phase transition, is more 
than adequate to demonstrate the underlying mathematics of the BGP. 

4.1 Homogeneous Variant 

The homogeneous version of the model has farms in all cells. An initial, 
randomly chosen cell is infected and begins the simulation. We chose to 
mitigate edge effects by inserting the initial infection within (100,250) instead 
of using periodic boundary conditions, because the former is computationally 
far simpler while achieving the same end result of avoiding broadening of the 
phase transition. In each run of this homogeneous variant, the transmission 
probability per day 11 is constant. Recall that 11 is applied only to healthy 
farms. 500 runs were averaged for each value of 11. Results were also obtained 
using 1000 runs and showed no significant differences from the 500 batches. 
The model parameters of time-lag (number of days between end of incubation 
and beginning of culling) and incubation were set to two days and one day 
respectively. The real FMD would have these as variable. 

In Fig. 5, the ratio of the number of surviving farms (not culled) to the 
number (homogeneous case: 90,000) initially present is plotted as a function 
of n. A phase transition is indicated. It represents the percolation critical 
point whereby the disease has either zero density or approaches 100% density 
(as one goes above the critical probability He) in the limit of very large 
landscape area (limit L oo). In the two-dimensional cellular automata 
lattice used here, above and very near the critical probability He, the disease 
density will scale as (11 — YlcY^^^; as soon as the conditional probability 
moves away from this scaling region, the disease density 100%. The 
critical transmission probability He for the homogeneous case is ~ .33. 
The epidemic time durations versus 11 are plotted in Fig. 6. Again, a peak 
is seen at 11 Ri .33. In order to be certain this is a phase transition, we run 
the homogeneous variant using a smaller size lattice of 100 x 100 cells and a 
larger size lattice of 700 x 700 to compare directly to the 300 x 300 case. Fig. 
7. As expected, the smaller lattice has a pronounced widening and the larger 
lattice has a pronounced narrowing. We thus have an unmistakable phase 
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transistion. In terms of the 300 x 300 lattice used, this critical probability is 
reached for an average cattle density ~ 1500; since pigs have a much higher 
FMD susceptibility coefficient, this critical probability is achieved with a 
much smaller number of pig population, or with a smaller combined pig and 
cattle population. 

4.2 Heterogeneous Variant 

The heterogeneous model differs from the previous homogeneous case by 
randomly selecting some cells to be vacant. As mentioned earlier, vacant 
cells are an impediment to disease diffusion. Runs were done on the same 
300x300 lattice and the farm occupancy factor, f, was chosen tohev — 0.75. 
In Fig. 5 we plot the survival ratio as a function of 11, while in Fig. 6, the 
epidemic duration as a ftinction of 11 is plotted. Again, a clear signal of a 
phase transition is present with the critical conditional probability per day 
transmission of 11^ ~ 0.44 for this heterogeneous case. 



4.3 Analytical Solutions for Homogeneous and Hetero- 
geneous Variants 

The homogeneous and heterogeneous models are amenable to computation of 
their critical conditional probabilities per day. Let us recalP that the square 
lattice has the critical (site) percolation density of 

Ppercolation — 0.59273 (2) 

For our situation here, the critical percolation density enters into the prob- 
ability of finding animals in a cell, allowing the continued diffusion of the 
infection. Cellular automata models have different neighborhood rules which 
cause the critical value of parameters to vary from model to model. In the 
case here, the diffusion probability P is a function of Ppercolation- We ask, 
"What is the probabihty of the disease spreading out from IP Cell[i][i] when 
nearest neighbor and next-to-nearest neighbors are culled?" 

Refering to Fig. 2, all farms with diffusion probability P or P/3 are 
culled, so the disease starting out from IP Cell[i]|j] must reach cells further 



11 



■ homogeneous 
A heterogeneous 





I 


■ A i 
■ 


k 

A 












■ 


i 












1 


1 














I 


▲ 
▲ 












1 

V 

1 


I.' 


k 

1 





0.1 0.2 0.3 0.4 0.5 0.6 0.7 
Probability per day 



Figure 5: Farm survival ratio (fraction of farms not culled) versus conditional 
transmission probability per day, 11, for both homogeneous and heteroge- 
neous cases. 



out within the alloted time-lag. The time-lag of 2 used here means that 
if IP farm A has infected a CP or DC, designated by farm B, on the first 
day of contagiousness, the disease has one additional day to use farm B as a 
new source of infection. We first compute the diffusion probability that the 
disease can go outside the culled area. 

In the first day of contagiousness, the only cells reachable outside the 
future culled those that have diffusion probability of P/9 and there 

are 8 of them, so first day diffusion P first day is 

P - ^ 

J- ftrstday g ■ y'-'J 

However, if the infection successfully infects a CP or DC on day one, it has a 
much better chance of escaping outside the culled area on day two. On day 
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Figure 6: Epidemic duration time (in days) versus transmission conditional 
probability per day, 11, for the homogeneous and heterogeneous cases, show- 
ing a peak duration which indicates a phase transition at a critical transmis- 
sion probability. 

two, each of the first day's successful jumps is itself a new Cell[i][j], which 
spreads out according to Fig. 2. Thus we need to keep track of the second 
day's diffusion relative to the future culled cells. So, for example, if on day 
one, it reaches a cell with probability P, and on the second day reaches its 
neighboring cell of probability P, the disease only reaches a future culled 
cell, so this mode adds nothing. To achieve diffusion, the second day jump 
from farm B must reach at least its P/3 cell. The computation of all possible 
routes is simple, but tedious. In Table 2, we tabulate the probabilities. The 
last column is computed by multiplying the first day's probability by the 
second day's probability and summing over all of the routes that avoid the 
future culled cells. 



13 



tn 
>> 

(0 

2, 

c 
o 



900 
800 
700 
600 
500 
400 
300 
200 
100 






■ 


700 lattice 








• 


100 lattice 








▲ 


300 lattice 
























■ 




























■ 














1 


1 










■ 




■■ 


■ 








■ — 












1 






i 




i 






i 














^ 








• 





0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 
Conditional probability per day of infection 



Figure 7: Homogeneous epidemic duration time (in days) versus transmission 
probability per day, 11, for the 100 x 100 lattice, the 700 x 700 lattice and 
the original 300 x 300 lattice 



By adding up all the probabilities and seting the sum to equal one, we 
define the meaning of the critical diffusion (percolation) probability Pc 

9 9 ^ ^ 

Had we used a time-lag of three days, then the resulting polynomial would 
have been cubic, changing to quartic for a time-lag of four days and so on. 

The infection cannot continue to percolate, however, unless the cell it 
has jumped to actually has animals. These potential jumped cells may have 
already been culled in earlier time periods. Futhermore, a likely jumped cell 
may be originally vacant, with probability 100 x (1 — z/), as happens in the 
heterogeneous case. The probability Panimais that any of the possible jumped 
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first probability 


second probability 


total probability 


P/9 
P 


P 


8P/9 



P 
P 

P/3 
P/3 
P/3 


P/3 
P/9 

P 
P/3 
P/9 
P/9 


20PV3 
16PV9 

8PV3 
32P79 
48PV27 

8P/9 



Table 2: Basic probabilities to avoid culled cells 



cells has animals at exactly the critical diffusion limit is 

Panimals — ^Ppercolation (5) 

Thus the critical diffusion probabihty Pc for disease percolation is a prod- 
uct of the critical transmission probability per day for disease transmission, 
lie, when animals are present, times the probability that animals are actually 
present: 

Pc X Z/ X Ppercolation ■ (6) 

Solving Eq.(3) gives 

He X i/ X Ppercolation = 0.1984 . (7) 

Letting (Tmuitipie — 1/.1984, where each infection generates a possible (Tmuitipie 
new ones outside of the culled area, the fundamental equation for disease 
percolation becomes 

lie X CTrnultiple X Z/ X Ppercolation — 1 (8) 

The product He x z/ x Ppercolation is the probability that there are animals 
in a cell and that they will obtain the disease. The number of cells outside of 
the future culled area succesfuUy jumped to is (Tmuitipie- Their product must 
be equal to one in order for the disease to be sustainable. Those readers 
familiar with the design of nuclear reactors will see that Eq.(8) is the neutron 
equation for criticality. 
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The analytical solution of He for the homogeneous and heterogeneous 
cases is given in Table 3. For the heterogeneous runs, the constant entering 
the landscape farm density probability, u, had the value of 0.75. 

Eq. (8) explains why there are two phases for the BGP: if the average 
diffusion probability is less than that required for escaping the culled cells, 
then the infection cannot be sustained, and damps out; in the limit of large 
landscape, the disease density is zero. Unfortunately with FMD, the long 
spatial tail and time-lag require an enormously impractical culling area. 



homogeneous case 


heterogeneous case 


0.335 


0.335/z/ 



Table 3: Analytical critical conditional probabilites per day He 

In variant 3, where a stochastic conditional probability per day is used, 11 
is required to satisfy the initial condition Hlow < H < Umax, for each farm. 
That is, each Cell[i][j] that happens to contain a farm (probability 100 x u 
that a cell contains a farm) has assigned to it a random 11, which is the same 
as saying that the Cell[i] [j] has a random animal population. 

4.4 BGP Variant 

This variation of the model is the closest to the BGP. For each run, the 

location of farms on the 300 x 300 lattice is random, as in the previous 
heterogeneous case (again the probability that a cell contains a farm is 100 x 
u). In addition, the animal population of each farm varies randomly, but 
fixed for that entire run. Thus each run represents an idealized realization of 
England with an initial random location of farms, each with an initial random 
farm animal population, categorized by a cell's value of 11. This 11 initial 
information is stored in an extra array: Prob[i] [j] holds the randomly assigned 
n for Cell[i][j], if this cell has a farm. When an IP Cell[i]|j] is identified, the 
location of this IP to other cell locations Cell[i'] [j'] containing healthy animals 
is determined, and the probabihty Prob[i'][j'] is then used (refering back to 
Fig. 4, this is the step "if true, roll probability ..."). Thus, for example, if 
Cell[i][j] is an IP, then Prob[i-l][j+l]/3 is used for the probability that a farm 
containing healthy animals in Cell[i-l][j-|-l]becomes infected. 
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What is there left to vary? FMD has the frightening abihty of non-local 
disease transport, aka the lightning factor. As long as there is one IP on the 
landscape, there is a lightning factor for every cell containing healthy animals. 
This lightning factor is a constant probability per day of non-local infection 
of healthy animals, a purely random effect. Examples of the lightning factor 
for FMD include airborne phimc of virus, and the spread of FMD to locations 
outside the local area by movement of persons, vehicles and animals through 
markets and agents. In Table 3, we present the results. The lower value 
of n (0.05) represents about 500 cattle (0.05 probability for 1 mile distance 
from an IP), while the higher value of 0.70 represents about 2500 cattle (0.70 
probability for 1 mile distance from an IP). 



Plow 


Pmax 


Survival ratio 


Duration (days) 


Lightning 


Farm density 


0.05 


0.70 


0.932 


134 





0.75 


0.05 


0.70 


.122 


386 


0.0001 


0.75 



Table 4: Variant 3 Results 



For zero lightning factor, the disease is in the low mortality phase. A 
lightning factor as low as 0.0001, however, pushes the disease into the high 
mortality phase. No attempt was made to analytically predict variant 3. 

5 Conclusion 

The underlying mathematical model of the 2001 BGP in the FMD epidemic 
possesses a phase transition. For a small average disease transmission kernel, 
the infection is quickly eliminated at hotspots and in the limit of a large 
area, the density of the disease is zero. A second outcome occurs if the 
average disease transmission kernel and its associated lightning factor are 
greater than well-defined numerical critical values. This results in the disease 
propagating throughout the landscape and, in the limit of a large area, the 
disease density approaches one as the conditional probability is raised above 
the critical value. FMD has wreaked havoc. The existence of the phase 
transition comes from the epidemiological feature of FMD: time-lags exist 
for any culling policy relying on the appearance of clinical symptoms. 
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In the United States, there are major cattle producing g. Florida) 

which have similar cattle density and spatial farm distribution as England. 
Thus implementing the BGP would be a great risk to the United States. The 
authors suggest the following alternative strategy: it is possible to detect 
the presence of the virus in animals before clinical symptoms or antibodies 
appear. This is done by PGR (polymerase chain reaction: a method to du- 
plicate RNA or DNA and thus amplify the amount of genetic material in 
the sample) amphfication of viral RNA in the host's serum. Recently, it has 
been shown^ that the FMD virus can be detected this way in animals. This 
technology means that the time-lag can be eliminated if mobile PGR labora- 
tory units can be made to return results in a few hours. Stationing of these 
units at ports-of-entry and within the interior will allow rapid first-responder 
capability. Identifying the virus at an IP before the clinical symptoms ap- 
pear means the culling policy will eradicate the disease efficiently: with no 
time-lag, no landscape percolation exists. Goincident with the PGR mobile 
units is the critical requirement of creating a FMD predictor model, using 
a complete database of farm animals. The FMD predictor model would 
direct the movement of the mobile units to predicted high risk localities. 
Without such a code, one squanders away the only real assest that will be 
effective in confronting this most formidable microorganism. With the use 
of first-responder units, a draconian movement restriction policy is essential 
to lowering the lightning factor of non-local disease transport. Finally, the 
United States must design and manufacture, as quickly as possible, mobile 
crematoria with high throughput. 
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